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ABSTRACT 



Several numerical methods used in the calculation of hydrodynami 
shocks were investigated. Particular attention was given to the 
artificial viscosity approach of Von Neumann and Richtmyer and its 
aoolication to the "PUFF" numerical scheme. The particle model 
approach of Ludford, Polachek and Seeger, and the method of Lax were 
also considered. 
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I. INTRODUCTION 



Frequently the numerical calculation of comoressible fluid flow is 
complicated by the presence of shock waves. The difficulties arise 
from the fact that shock waves propagate discontinuities in velocity, 
pressure, and other variables characterizing the fluid flow. Various 
approaches to this problem have been offered, each having its own 
desirable and undesirable characteristics. These approaches generally 
fall into two categories. The first category is related to the study 
of the viscosity of the fluid, whereas the second category is dependent 
on the conservation form of the hydrodynamic equations. The first 
category of approaches will receive primary attention. 
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II. NUMERICAL CALCULATION OF HYDRODYNAMIC SHOCKS 



USING AN ARTIFICIAL VISCOSITY FACTOR 
A. INTRODUCTION 

The equations describing perfect compressible fluid flow, in the 
presence of shock waves, produce solutions with discontinuities. 
Investigation of the physical situation shows that the true discontinu- 
ities, however, cannot occur due to the viscosity, or inner friction, 
of the fluid. These equations ignore the viscosity of the fluid and 
do not accurately represent the physical system. The addition of vis- 
cosity terms to this system of equations shows that the fluid behavior 
inside the shock region is nonlinear but continuous. The viscosity 
of the fluid is negligible outside of the shock and significant inside 
the shock. The original intent then, was to replace the shock reqion 
by a discontinuity and treat the shock as a two-sided boundary. The 
size of the jump discontinuity, or rather the boundary values, would 
be prescribed by the Rankine-Hugoniot equations. However this approach 
has several drawbacks. First, the presence of a discontinuity comoli- 
cates the use of a numerical scheme to solve the problem. Secondly, 
the shock wave is in motion, and hence, the boundary is moving. Thirdly, 
since irreversible thermodynamic changes of state take Diace across a 
shock region, an increase in the specific entropy of the fluid must be 
added to the original jump conditions. And, lastly, this approach does 
not represent the physical situation in that there is no indication of 
the behavior of the fluid inside the shock region. 

Since the addition of the true viscosity terms severely complicates 
the system of differential equations. Von Neumann and Richtmyer [Ref. 7] 
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have suggested the addition of a "pseudoviscosity" term to the equa- 
tions of non-viscous flow, which will not complicate the system as much 
as the true term would. This term must, of course, conform to several 
restrictions . 

B. THE BASIC EQUATIONS 

The equations describing one-dimensional flow of a compressible 
fluid are as follows: 



VU,t) ' 3X 


C2.1) 


u(x,t) = fi 


(2.2) 


aU _ 3 Cp+q ) 

P 0 3t 3t 


(2.3) 


- 0 


(2.4) 


3V _ au 

P 0 3t 3X 


(2.5) 


E = P V . , 

L wry 


(2.6) 



where x is the Lagrangean coordinate, X = X(x,t) is the Eulerian 
coordinate (i.e., X(x,t) gives the position, at time t, of a fluid 
element initially at position x), p q (x) is the initial density, V is 
the specific volume, U is the fluid velocity, p is the static fluid 
pressure, E is the internal energy Der unit mass, y is the ratio of 
specific heats (i.e., Cp / c y ), and q is the artificial viscosity. 
Equations (2.3), (2.4), and (2.5) are the equations of motion, of 
energy, and of continuity respectively. Equation (2.6) is the equation 
of state for a perfect gas. 
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C. THE ARTIFICIAL VISCOSITY FACTOR 



The expression for q must satisfy the following requirements: 

1. Equations (2.3), (2.4), and (2.5) must possess solutions 
without discontinuities. 

2. The thickness of the shock must be everywhere of the same 
order as Ax (the length increment) used in the numerical 
calculation. 

3. The effect of q on eqns. (2.3) and (2.4) must be negligible 
outside the shock region. 

4. As ax+0, the solution must approach a state with a jump 
discontinuity prescribed by the Rankine-Hugoniot equations. 

Apparently, these requirements are not enough to uniquely define q. 

In any case the expression that Von Neumann and Richtmyer developed is 



q 



(p 0 cax ) 2 

V 



aV aV 
at at 



(2.7) 



where c is a dimensionless constant near unity. By the use of equation 
(2.5) q can be written as 



q = - (cax)^ all . 

„ ax 



all 



ax 



( 2 . 8 ) 



Von Neumann and Richtmyer have proven that q satisfies the above 
requirements fora particular case of steady-state plane shock. They 
conjecture, however, that the artificial viscosity approach would be 
equally suited to more complicated multi -dimensional flows. The problem 
they consider is the example of a one-dimensional shock wave separating 
two regions of constant state. This simulates the situation that occurs 
when a piston is pushed at a constant velocity into a long tube con- 
taining a fluid initially at rest. After the shock has traveled a 
sufficient distance from the initiating piston, it moves at a constant 
speed, s. In the absence of an artificial viscosity term, the specific 
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volume, V, at some time t, is given by figure 1 Since we are con- 
sidering steady-state solutions only, the solutions depend only on a 



linear combination of x and t given by 

w = X - St. 

Define 


(2.9) 


M = p 0 s. 

Now 

U(x,t) -*• U(w) = U(x - st) 
implies that equation (2.3) becomes 


(2.10) 


M dU = d(p + q) 
dw dw 

since 


(2.11) 



3U _ 3U _ aw 3U = 3U 
aw p o s 3w p o at aw p o at 



3 ( p+q ) _ d(p+q) aw 

ax dw ax 

Similarly, equations (2.4) and (2.5) become 


_ d(D+q) 
dw 


3w + (p+q) Bw = 0 


(2.12) 


and 

m dV dU 

' M ^7 = ^ 

Then equations (2.11) and (2.13) qive 


(2.13) 


,,2 dV _ d(p+q) 
dw dw 

and equations (2.12) and (2.14) give 


(2.14) 


dE + CfpigM + m 2 dV . 
dw dw dw 


(2.15) 
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Von Neumann and Richtmyer then integrated equations (2.13), (2.14), 
and (2.15) with respect to w giving 



MV + U = C ] 

M 2 V + p+q = C 2 

E + (p+q) V + 1/2 M 2 V 2 = C 3 

as solutions of equations (2.13), (2.14), and (2.15) where C-j , 
and are constants of integration. Let the initial and final 
be denoted by: 

As w-*»; V-*V . , p->-p.j , E^E .. , q->-0 
As vi-*‘~ c °y V^V^, p y P^r , E+E _p , q^O 

Since V.. and V f are particular values of V, and p^. , and P f are 
values of p, they satisfy equation (2.17) giving 

M 2 V. + p. + 0 = c 2 

M 2 V f + p f + 0 = C 2 , 

which implies that 

M 2 (V.-V f ) = (p f -p.) . 

A similar argument using equation (2.18) yields 

E f + P f V f + 1/2 M 2 V f 2 = C 3 

- Ej - Pj V, - 1/2 H 2 V 1 . 2 = C 3 

from which it follows that 

(E f -E.) = 1/2 M 2 (V. 2 -V f 2 ) + p.V. - p f V f . 



(2.16) 

(2.17) 

(2.18) 

C 2 , 

values 

(2.19) 

( 2 . 20 ) 
particular 



( 2 . 21 ) 
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Then by equation (2.21) 





1/2 (p f -p 1 )(V i +V f ) + p.V 1 - p f V f 



1/2 (p.V. + p f V t - p f V f - p f V f ) , 



or 



(E f -E.) = 1/2 (p.+p f )(V r V f ) 



( 2 . 22 ) 



Von Neumann and Richtmyer point out that equations (2.21) and (2.22) 
are independent of q, providing q->0 as w-*±», and in fact are the Hugoniot 
equations. Requirement (4) is then satisfied since it has just been 
shown that the Rankine-Hugoniot conditions are satisfied for the flow 
sufficiently far from the shock region. Requirement (4) may be 
examined in an alternate fashion. Let Z(x,t) be the solution to the 
set of equations describing non-viscous flow. I.e., Z(x,t) has a jump 
discontinuity, prescribed by the Rankine-Hugoniot equations, in the shock 
region. Now let w(x,t,Ax) be the solution of the equations of viscous 
flow corresponding to a fixed q, or more accurately, a fixed ax. Then 
we require 



By the very fashion in which q was introduced equation (2.23) is satis- 
fied. Note that q actually has the dimensions of pressure and enters 



approach zero. Consequently, the system of equations describing viscous 
flow approaches the system describing nonviscous flow as our mesh size 
approaches zero. 



limit w(x,t,Ax) = w(x,t,0) = Z(x,t) (2.23) 

AX->0 



equation (2.3) and (2.4) in the forms and (p+q) respectively. 

Since q is continuous, |^>0 as q-*0, and hence, all viscosity terms 

o X 
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The question naturally arises as to why Von Neumann and Richtmyer 
have created a process whereby decreasing the mesh size causes w(x,t,Ax) 



to approach a solution that does not represent the physical system. 
The answer is a heuristic one in that it must be possible to make q 
arbitrarily small to accomodate arbitrarily thin shock waves. It 
should also be noted that q is physically artificial and was created 
for numerical convenience only. 

To investigate the shape of the shock Von Neumann and Richtmyer 
consider solutions satisfying 



Equations (2.24) are normally the situation characterizing a shock 
moving to the right. Then equation (2.7) yields 



nV aV 

g^iO, or equivalently, —> 0. 



(2.24) 



qV = (Mcax ) 2 {jjt) 



2 



(2.25) 



Now from equation (2.18) 



E + 1/2 M 2 V 2 = C 3 - ( p+q ) V 



E - 1/2 M 2 V 2 = c 3 - ( p+q ) V - mV 

= C 3 - V[ (p+q) - M 2 V]. 



And then equation (2.17) gives 



E - 1/2 M 2 V 2 = c 3 - c 2 v . 



(2.26) 



Then by equation (2.6) 



pV = (E)( y -1) 

= 1^1 mV + C 3 ( y -1) - VC 2 ( y -1) 
= [ ^ ]M 2 V 2 + c 4 v + c 5 



(2.27) 
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where 



-c 2 (y-1) 



(2.28) 



C 



4 



C 



5 



C3 (y-1) 



(2.29) 



Then equation (2.17) yields 

qV = C 2 V - pV - M 2 V 2 

= c 2 v - c 4 v - c 5 - [ ^ M 2 V 2 + M 2 V 2 ] 

qV = C 2 V - [^^-] M 2 V 2 - C 4 V - C 5 . (2.30) 

Now for V = V.. and V = V^, q = 0. Therefore since the right side of 
equation (2.30) is quadratic in V and vanishes for V = V. and V = V^, 

qV = M 2 V 2 + (C 2 -C 4 )V-C 5 

= M 2 (-V+V.)(V-V f ) . (2.31) 

Then equation (2.25) yields 

(Hoax) 2 [£] " ^ (V 1 -V)(V-V f ) 

<Cix) Z [£] 2 = (V 1 -V)(V-V f ) (2.32) 

To solve equation (2.32), Von Neumann and Richtmyer proceed as follows: 
Let, 

V,+V f V.-V f A 

A = V - , A 0 = -4" 1 > B = £ (2.33) 

Then, 

CAX £ = ^ V2 ( Y v ) 1/2 ( v - v f ) 1/2 
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Therefore, 



or 



giving 



Hence, 



where 



Finally, 



- A]’/ 2 [A - 

- C(a 0 -a)(a+a 0 )] 1/2 . 



C4X 35 * Pj4 172 [(A 0 2 -A 2 )]' /2 , 



(2.34) 



A 2 -A 2 



c4x( r-* 35 = Pr] 1/2 t- 2 -^ 172 > < 2 - 35 > 

0 



C4X £ = P^ 1 ] 172 [1 - b 2 ] 1/2 . 



(2.36) 



dp 

dw = 72 “x t 



w ■ [r|r] ,/2 c4x/ 



dB 



w = 



arcsin B , 
o 



(2.37) 



-0 ' ^] 1/2 cax 



(2.38) 



A = A B = A sin Jf 
0 0 w. 



V.-V- 

— — - [sin jj-], (2.39) 
2 w o 
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or 



V V f w V-+V- 

V = — — - sin J$- + — — — . (2.40) 

2 w o 2 

V is obviously continuous and hence requirement (1) is satisfied. 

Von Neumann and Richtmyer state that w Q is a measure of shock thick- 
ness. Thus, if c is near unity, w Q is 0(ax) and requirement (2) is 
satisfied. 

Now taking the derivative of V with respect to w and setting it 
equal to zero gives 



dV 

clw 



V v f 1 

2 v/ 



cos 



w_ 

w„ 



0 . 



(2.41) 



Hence 



cos 



w_ 

w. 



= 0 , 



(2.42) 



_ (2n+l )tt s (n an integer), as points where V assumes 



giving w ^ - 0 

its relative maxima and minima. Similarly, 



d 2 V , 

TT =( 

dw 



• V v f . i 

9 ) 



W 



i - 0 

w o 0 



(2.43) 



impl ies 

sin = 0 , (2.44) 

0 

giving w = nirw Q , (n an integer), as inflection points for V. Now 

for w = - 7 w , 

2 0 

V V V V f 

v = V V f Sin (-f) + = V f . (2.45) 

2 2 
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V 



sin ( 75- ) + 



(2.46) 



For w = vj , 

c. 0 



V v f 



v , tv f 



v t 



And for w = 0 

V i +V f 

V = - L — 1 - (2.47) 

2 

Finally, since only solutions satisfying ^ 0 were considered, V is 

non- oscillatory. As a result of equations (2.41) through (2.47) 
figure 2 represents our solution. 

Note that for this particular problem of steady-state plane shock, 
q=0 outside C-f w 0 > ^w 0 ](i.e., the shock layer) since = 0 in this 

region. Normally, outside the shock region, q would be negligible 

p 

in comparison with the static pressure P because of the factor (ax) 
in equation (2.7) and a relatively small value for ^r. However, inside 

the shock layer q is comparable to p because of the abnormally large 
value of encountered in that region. Hence requirement (3) is 

satisfied. Therefore, for this particular case of steady-state plane 
shock, Von Neumann and Richtmyer have shown that their expression for 
q conforms to all the necessary restrictions. 



D. STABILITY OF THE DIFFERENTIAL EQUATIONS 

Our next concern will be the effect the introduction of an artificial 
viscosity term has on the entire system of differential equations. 

Before investigating the stability it should be noted that much of 
the computational work actually done will be omitted for the sake of 
brevity. Instead, reference will be made to the approach and ideas 
involved. 
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On a given solution U(x,t), V(x,t), etc., a small perturbation 
6l), 6V, etc., is superimposed. Then the system will be stable if the 
perturbation can be kept arbitrarily small for all t >_ T Q by initially 
choosing the perturbation at time T q to be sufficiently small. There- 
fore consider the following variations: 

U -*• U + 6U 
V + V + 6V 
p -*• p + 6p 

q -> q + 6q 

We then obtain our equations of first variation (i.e., higher order 
variations are considered to be negligible): 

P n 3(6U) _ 3(6p+6q) 

°~1T ‘ " 3X (2.48) 



||tY4p+(y-l)«q] + [ T p+(t-l)q] + V 



dt 



at 



+ H <V - 0 



(2.49) 



, n _ (cAx) 4 " 9 U 

6q - o — — 

T 



3X 



3U 



3x 



6V - 2 



(CAX)‘ 



9U 



9X 



9(aU) 

9X 



(2.50) 



9(fiV) _ 9 ( 6U ) 



9t 



9X 



(2.51) 



Equations (2.48) through (2.51) are a set of simultaneous, linear 
differential equations in 6U, <sV, 6p, and sq. Their coefficients are 
composed of terms depending on the solution functions U, V, p, and q. 
Since U, V, p, and q are considered to be smooth, well-behaved functions 
of x and t, they will be considered to be constants in a small region. 
Equations (2.48) through (2.51) were combined into one equation and a 
separation of variable technique was employed. Using this technique 
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possible solutions were of the form: 



(2.52) 



6 U = <5U o e at (cos kx + i sin kx) = 6U 0 e 1kx+at 

6V = <sV o e at (cos kx + i sin kx) = 6V 0 e 1 ' kx+0 ‘ t 

6p = <SP o e at (cos kx + i sin kx) = 6p o e ll<x+at 

<5q = <$q 0 e at (cos kx + i sin kx) = 6q Q e ll<x+at 

where 6U . 6V . 6p„, <5q . k, and a are constants and k is real. Sub- 
0 0 0 0 

stitution of equations (2.52) into equations (2.48) through (2.51) 
yields the following set of simultaneous homogeneous linear equations 



in 6U q , 6V q , 6p 0 , and 6q Q : 

RS = 0 

where R is the following 4x4 matrix: 

Ik 



[| V 



(2.53) 




0 

0 




(2.54) 



r6l) e 
0 



i kx+at\ 



ikx+at 



ikx+at 



(2.55) 



.ikx+at 



20 



Now if equation (2.53) is to hold for a nontrivial S, then 



DET(R) = 0 

The characteristic equation is: 



(2.56) 



, n 2 9V . p o a 3V av v 2 
(ap o ) Y 3t 2 ~Jt (kcAx) 



9U 

3X 



' 3 If (kc4x) ‘ 



3U 


3U 


3X 


3X 



+ k a[yp+(Y-l )q] 



+ p + 2p^a^(kcAx)^ 



3U 



3X 



V <k“x) 2 g 



3U 



3X 



If ■ ° 



(2.57) 



Equation (2.57) establishes the relationship between a and k. By fixing 
k and examining the corresponding a, we can investigate the behavior of 
the perturbation. It should be noted that the behavior of the pertur- 
bation must be examined both in the shock regions and in the normal 
regions. In the shock regions all terms will be retained. In the 
normal regions, terms containing dissipative factors (i.e., terms con- 
taining ax) are dropped. Now we are interested only in the a*s corres- 
ponding to very large k's. Hence we retain only the dominant terms, 

in a and k, of equation (2.57). The dominant term in a is p Q a V. 

2 

The dominant term in k contains k and is therefore dominated by 
3U 



2p Q a^(kCAX)^ 



3X 



, which is the dominant term in ak. Hence in the 



shock regions, equation (2.57) reduces to: 



p **a^V + 2p a^(kCAX)^ 
0 0 



3U 



3X 



= 0 , 



(2.58) 



giving 



a = 



-2(kcAx)‘ 



3U 

3X 



(2.59) 
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and in the normal regions 



k 2 ayp 



2 3., 
+ p Q a V 



= 0 , 



(2.60) 



giving 



2 

a 




(2.61) 



Thus, in the normal regions our original system of differential 
equations is stable, and in the shock regions, the system is asymptotically 
stable. Von Neumann and Richtmyer also point out the terms in the 
equations of variation that lead to the dominant terms in equation 
(2.57). In the shock regions: 



3 ( 6 1) ) ~ _3 2 (sU) 

~ at" ~ 0 „ 2 

3X 



where 



a 



2(cax) 2 

V ^o 



3U_ 

3X * 



and in the normal regions: 



(2.62) 



(2.63) 



3 2 (5U) ~ c 2 3 2 (6U) 
at 2 0 3X 2 



(2.64) 



where 



S 



2 

o 




(2.65) 



E. FINITE DIFFERENCE SCHEME 

To solve the system of differential equations several finite 
difference schemes could be used. The one Von Neumann and Richtmyer 
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offer is ingeniously simple. The central differences used are skill- 
fully staggered, taking advantage of the artificial nature of q. 

Let our rectangular mesh with increments ax and At, and integers 
i =0, 1, 2, ..., L; n=0, 1, 2, ... be contrived in the following 
fashion: 



U i+l/2 = u [(i + l/2)Ax, nAt] 



V' 



i+1/2 



= V[(i+1/2)AX, nAt] , etc. 



(2.67) 



After rewriting equation (2.4) the finite difference scheme is as 
follows: 



p o 1 



Un+1/2 _ ^n-1/2 



_n + r, n_1 /2 _n 
p i+l/2 q i+1/2 ~ p i-l/2 



At 



AX 



- d "- 1 / 2 

q i-l/2 



( 2 . 68 ) 



AX 



V n+ 1 _ v n 

V i+l/2 v i+l/2 
At 



.,n+l/2 , .n+1/2 

U i+1 ~ u t 

AX 



(2.69) 



q 



n+1/2 

i+1/2 



2(cax) 2 (U 



n+1/2 

1+1 




(ax) 



2 



(v?+ 



1/2 



+ V 



n+1 

i+1/2 



) (2.70) 



^n+1 , n , i tr\ \/fl \ 

[ Y P i+1/2 P i+l/2 + (Y-DqJ+]/2 ] (V i + V^ " V i+l/2 } 



At 



+ 



A ,n+1 , .,n wJi+1 n \ 

(V i+l/2 V i+l/2 )(p i+l/2 ' p i+1/2 ) 

2At 



0 (2.71) 



Since central differences were used the discretization error will 
be 0(ax) 2 and 0(At) 2 . 
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F. STABILITY OF THE FINITE DIFFERENCE SCHEME 

Having shown the original system of differential equations is 
stable and having chosen a suitable finite difference scheme one must 
now show that the difference equations are in fact stable. A finite 
difference scheme is said to be unstable if the rounding error, intro- 
duced in approximating the numerical solution, grows exponentially with 
each iteration, making nonsense of our numerical data. This concept is 
analogous to the one for differential equations in that the rounding 
error corresponds to a small perturbation in the numerical solution. 

Von Neumann and Richtmyer have shown equations (2.68) through 
(2.71) to be conditionally stable (i.e., stable only for certain com- 
binations of ax and At). Hence certain restrictions will be placed 
on the choice of ax and At. It should be further noted that these 
restrictions will not be the same in shock regions as in normal re- 
gions. As before much of the computational work actually done will be 
omitted for brevity and clarity. 

Outside shock regions the stability criteria of Von Neumann and 
Richtmyer is 

A+ n 1/2 

£(1V) <1. (2-72) 

e 0 V 

Equation(2.72) is the usual stability criteria encountered when hydro- 
dynamic equations of the form of equation (2.64) are approximated by 
central differences. 

A similar analysis in the shock regions yields 

< 1. (2.73) 

(ax ) 2 “ 
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From equations (2.38), (2.39), and (2.62) 



= scax [ ■ 1 y f ] [lii] 1/2 cos w 



w 



Equation (2.73) then becomes 



J + 1 4- w 

.. ? 1/2 J-l + sin w 

— <( _2) o 

ax - v y+r 



. w 

4 sc cos — 

w o 



where 



Now let 



V, 



J = 



f O 

o 



J+1 . w 
O^T + sin w" 



cos 



w_ 

w. 



for -tt/2 w Q < w < ir/2 w Q . 



J+1 w . . w 

TT sec — + tan y- , 

J - ‘ w ~ 

0 0 



Then, 



F' - sec tan £- + sec 2 

J- I VT W. W_ 

0 0 0 



W 

J+1 5in w 0 ,_J 

^ cos 2 £- cos 2 £- 

0 0 



2 w r / J+1 \ a w , *i*i 
* sec irl ; (xT)^'' i r+ 1] 
0 0 



Now setting F 1 = 0 implies 



(2.74) 



(2.75) 



(2.76) 



(2.77) 



(2.78) 
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(2.79) 



* W / J - 1 \ 

sin vT ' (j+T^ 



since sec e is never zero. Hence, we have obtained as a critical value 



J+1 J-l 
F = J-1 ~ J+1 

[i -^4 ] 1/2 
(j+i) 2 

ojl/2 

= (2.80) 

Noting that as •* \ from the left F(^-) -»• °°. Since only one critical 
o o 

point has been found, equation (2.80) must be the minimum value for 
F(~“) in the shock region (-ir/2 w < w < tt/2 w_). Consequently, 

W _ 0 0 

0 

equation (2.75) may be rewritten as 





1/2 



J 1/2 

^ 2sc(J-1) ^ 



From equation (2.21 ) 



M 



r'Pf-Pj 

[ w 



■] 



1/2 



Equations (2.22) and (2.6) give 



(2.81) 



(2.82) 



so that 



Pf v f . Pi v i = 1/2 ( Pl + Pf )(V r V f ) 

Y~1 yl 

= 1/2 p f (V.-V f ) + 1/2 p f (V.-V f ) , 

(2.83) 

1/2 p f (V i -V f ) = [1/2 + yyOp.jV.. 2 — » ( 2 - 84 ) 
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yielding 



5* - p f < v i- v f) ■ Ow p, • 

Y-l 

Then equations (2.10) and (2.21) give, after eliminating p.. , 

1/2 



(2.85) 



or 



S = 





p f ■ 




X 


1_ 1 


L<£X v f J 


V 


p o 


<VV 


s 



S = — s 



n rl 1 (xtlHlidi-j 

S ~ y (v+1)J-(y-D" J 



j-i 



1/2 



Finally 



[ 



J(y+1 ) - ( y- 1 ) 



-] 



1/2 



P-f 

Cy 



1/2 



1/2 



[ 



j( Y +l )-( y-1 ) 



] 



’of 



( 2 . 86 ) 



where S Q ^ is the speed of sound behind the shock, relative to x. 
Now equations (2.86) and (2.81) give 

S/2 



At < 1 



ax — 2S Qf c 



[ 0- 



(y-1 ) 
Iy+T7 



■]0 



(o-i V 



(2.87) 



Since V^, and hence 0, is generally unknown until the problem has been 
solved, the term inside the radical is replaced by its minimum with 
respect to J: 
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Let 



Q(J) 



[ (J - *yy ) J ] 
(J-D 



1/2 



( 2 . 88 ) 



Then 



q- = 



- 1/2 1/2 
1/2 (J-1 ) [( J-^yy ] • [ 2J-^|] - [( J-;^yy ] 



1/2(J-1 )(2J-^-) - [ J-^fy] J 



(J-1) 2 C(J-^|)J] 



(J-1)' 
J-: 
1/2 



(2.89) 



Setting Q ' = 0 gives 



0£1/2(£}>- 1] + l/2(^)- 0 



so that, 



But 



J = 



1 < J < 



nl 

y+3 

Y+1 

- 7T 



(2.90) 



(2.91) 



(2.92) 



since J = ^y corresponds to an infinitely strong shock and J < 1 
implies that Q is imaginary. Since yJj- < 1, there are no critical 
points in the interval [1, ^yy]. Hence J = ^yj is the point where 
Q assumes its minimum since Q can be made arbitrarily large as J-»-l 
from the right. Hence 



k1 h (i1 i )_ 

Y- I Y- I 



1/2 



'min 






Y 



1/2 



(2.93) 
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Finally substitution of (2.93) into (2.87) yields 



At < 

AX - 



1/2 



(2.94) 



2S of c 



Equation (2.94) is then a sufficient condition for stability of the 
difference equations inside the shock region. It should be noted that 
Von Neumann and Richtmyer have ignored boundary conditions in their 
stability analysis. Should derivative boundary conditions enter the 
problem an appropriate stability analysis will have to include a 
discussion of the difference equations used to approximate the boundary 
conditions. 
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III. INTEGRATION OF THE EQUATIONS OF TRUE VISCOUS FLOW 



An alternate method to the numerical calculation of hydrodynamic 
shocks is offered by Ludford, Polachek and Seeger [Ref. 5]. Their 
approach is similar to that of Von Neumann and Richtmyer in that their 
analysis examines the viscosity of the fluid. However, by replacing 
the fluid continuum by a particle model, they show that the true equa- 
tions of viscous flow can be numerically integrated. 



A. THE BASIC EQUATIONS 

The equations of one dimensional flow of a perfect viscous compres- 
sible fluid may be written 



where 



DU 
p Dt 


= iE. + 

3x 3x ’ 


(3.1) 


D£ 

Dt 


3U 
p 3x 


(3.2) 


„ 3U 
3x 


LO l+J 
0|0 

1 — 

Q_ 

II 


(3.3) 


P = 


PRT ; R = C V ( Y -1), 


(3.4) 


a = 


4 3U 

3^ 3X 5 


(3.5) 


S = 


C y log(p/p Y ). 


(3.6) 



T, a, S, and u are the temperature, viscous stress, entropy per unit 
mass, and coefficient of viscosity respectively. D represents the rate 
of change of the quantity written after it, if we move with the gas 
particle (Lagrangean viewpoint). All other quantities were defined 
previously. Now p may be eliminated from eqn. (3.3) by use of eqns. 
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(3.4), (3.6) and (3.2) in the following manner: 



a 



9U 

3x 



<*{| - 



RT 



DS 

Dt 



giving 



3U = p DS 
3X C v (y-1 ) Dt 



(3.7) 



Since 



we have 



Then 



Dt = C v Dt ^°9 (p/p Y )} 



{ 1 DR _ x Dfi) 
1 p Dt p Dt' 



(3.8) 



9U = J2_ / 1 D£ _ x D£ , 
3x y-1 1 p Dt p Dt 1 



(3.9) 



a 



9U 

8X 



1 

rpr 





(3.10) 



finally giving 



Dt = If [a ^ _1 ) ’ 



(3.11) 



Note that eqn. (3.11) is independent of p. 

Now Ludford et al . approximate the motion of the fluid by the motion 
of a system of particles. They consider flow inside a tube of unit 
cross-section. The tube is considered to contain (N+2) particles, each 
of mass m, and each moving along a fixed straight line. Each particle 
then represents that mass m of gas which initially had the particle at 
its center (see figure 3). Hence the forces encountered within the 
original gas may be approximated by the interaction of forces between 
the particles. 

j. 

If x^ is the x coordinate of the n u particle, and a dot represents 
differentiation with respect to time, the viscous interaction between 
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the n th and (n+1) th particles is given by 



'n+% 



_ 4 r X n+l - x n n 

" T P y j > 

J x n+l x n 



(3.12) 



as a result of eqn. (3.5). Now from eqn. (3.1) 



P x , 



"^ p n+?g " p n-%) ^ g n+% " a n-%^ h 

( V% - " X n-V> ’ * 



giving 



mx = “ (p n + % - p n-% } + - <W 



(3.14) 



since 



m = p( v% - W (1) - 



(3.15) 



The pressure interaction between the n th and (n+l) th particles is given 



by 



x n+l ' x n 



[(Y-Da n+ . - YP n .J. (3.16) 



n+% x n+1 - x n LV( ,/u n+% ,K n+% J 



If the constants u and y are those of the original gas equation then 
eqn. (3.2) is automatically satisfied. Hence eqns. (3.12) through 
(3.16) are sufficient to characterize the particle model. 

Before choosing a finite difference scheme, Ludford, Polachek, and 
Seeger nondimensionalize eqns. (3.12), (3.13), (3.14), and (3.16) by 
the following transformations : 

1 2 

p = pP = — p a 7 

H H o y M o o 




x = &X 

* 




(3.17) 
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qivinq in dace of the original equations of narticle motion, 



J n+% 



x ; + i 



- X n 



X n+1 X n 



X n" 



■ - (P n + % - V%> * A < E n + % ’ W 



(3.18) 



P' 



n+k 



X n+1 ' X n 



* n+1 ' - X n • 



where primes denote differentiation with respect to T, and 



*- 4 >^ [x(N+,)]1/2 



(3.19) 



B. FINITE DIFFERENCE EQUATIONS 

The finite difference scheme chosen to represent eqns. (3.18) is 
as follows: 



Y m+1 Y m Y m+1 Y m 

m+> 2 = _2 r- ( A n+1 ' Vi) - (*n ~ *n)-, 

h n+k ~ AT L / v m+% . v m v / v m+l , V nv 

U n+1 + Vl J " U n + n ' 



z 



m 

n+% 



1 rv m+ ^ 

2 L V% 



+ 




(3.20) 



„m+l 



2X m + X™- 1 
n n 



(4T) 2 [-(P™ +% 




) 



+ A(z 




-Z 



m 

n- 



& 



p m+l p ra 
v&\ n+% 



DKy-DAP" 1 ^ 



m+k 

n+% 



X 

2 



(P 1 



m 

r\+h 



P m+1 )l 

K n+V J 
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m th 

In the system, the quantity is the pressure between the n tn 

and (n+l)th particles at time T = mAT. It should be noted that the 
above system cannot be solved explicitly since the superscripts of T. 
are staggered. Hence eqns. (3.20) must be solved by an iterative pro- 
cedure at every AT/2 time increment. 

The stability analysis of the finite difference equations is quite 
similar to that of Von Neumann and Richtmyer. Equations (3.20) are 
numerically stable if, in the normal regions, 



AT < [ 





) 1/2 
-] 



(3.21) 



The system is unconditionally stable in shock regions. 

The aoproach taken by Ludford, Polachek, and Seeger has some 
advantages over that of Von Neumann and Richtmyer. First, in the 
particle model approach, the stability requirements are somewhat less 
stringent. Secondly, this approach gives a fairly accurate representation 
of the behavior of the fluid in the shock regions, whereas the artificial 
viscosity approach does not represent the physical system in the shock 
layers. However the approach by Ludford et al. does have a major dis- 
advantage. If the fluid under investigations has a very low viscosity 
an inconveniently fine mesh will have to be used in order that the finite 
difference approximations are accurate beyond the shock front. Otherwise 
AX will be larger than the thickness of the shock wave. It should be 
recalled that in the artificial viscosity approach q automatically 
adjusted to the shock thickness. 
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IV. THE METHOD OF LAX 



It was stated earlier that approaches to the numerical calculation 
of compressible fluid flow generally fall into two categories. The 
first category consisted of methods which examined the viscosity of the 
fluid. The motivation for this approach was that when a nonlinear term, 
multiplied by a small coefficient, is introduced into a differential 
equation, it may produce large changes in the behavior of the solution. 
The second category consists of methods which examine the basic equations 
of nonviscous flow, while allowing discontinuous solutions. An ingenious 
method developed by Peter Lax [Ref. 4] belongs to the second category. 

It is evident that the approaches discussed so far are heuristic in 
nature and lack theoretical preliminaries. Lax succeeded in makinq 
several theoretical observations which could tie together the various 
methods discussed in this paper. Unfortunately several of his observa- 
tions have been proven only for particular cases. 

Observing that all nonlinear systems of fluid dynamics satisfy 
certain "conservation laws", Lax states that any hydrodynamic system can 
generally be brought into the form 

U t + F x + B = 0, (4.1) 

where U is a column vector of unknown functions, F is a column vector 
such that F = F(x,t,U), and B is a vector coefficient. U is said to be 
a weak solution of eqn. (4.1) with initial value 4 if the integral 
relation 



/ / {lO + W F - WB) dxdt + / W(x,0)$(x)dx = 0 (4.2) 

0 L X 

-oo -oo 
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holds for every test vector W which has continuous first derivatives and 
vanishes outside of some bounded region in the x,t-plane. Equation (4.2) 
is obtained by multiplying eqn. (4.1) by W and then integrating by parts 
as follows: 



00 oo 



0 = / / {WU. + WF + WB} dxdt 

0 -oo t x 



/ [WU] dx - / / UW dxdt + / [WF] dt 

-oo 0 0 -00 0 -00 



OO CO 



/ / FW dxdt + / / WBdxdt 

0 -oo x 0 -oo 



(4.3) 



= h - l 2 + l 3 + l 5 • 



Now 13=0 since W vanishes outside some bounded region. Similarly 



I, = - / W(x,0)$(x)dx. 



(4.4) 



Hence we obtain eqn. (4.2). 

Clearly, the only requirement for U to be a weak solution is that 
it be continuous almost everywhere so that the (Riemann) integrals in 
eqn. (5.2) exist. Hence weak solutions need not be differentiable. 

The motivation for developing the notion of weak solutions is that in 
physical systems we are concerned with discontinuous functions that 
satisfy eqn. (4.1) almost everywhere. 

The Ranki ne-Hugoniot equations now have an interesting interpretation 
in terms of weak solutions; if U-j and are two genuine solutions of 
eqn. (4.1) whose domains in the x,t-plane are separated by a smooth curve, 
the two taken together will constitute a weak solution if and only if the 
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slope m of the separating curve and the value of U-j and on the curve 
satisfy 

1 (U-j - U 2 ) = F ( U-j ) - F(U 2 ) (4.5) 

Equation (4.5) corresponds precisely to eqns. (2.21) and (2.22). 

Havinq generalized the concept of a solution to a differential 
equation, we night expect some of the properties of the original concept 
to be lost. For instance, initial values do not in qeneral determine a 
unique weak solution. This property raises immediate difficulties. In 
nature nonviscous flow is described by a unique weak solution to the 
hydrodynamic equations, qiven an initial vector. Hence if our mathemati- 
cal model is a meaninoful one, there must be some other Drinciple that 
uniquely defines a weak solution. Lax offers some possibilities, the 
most likely beinq that the weak solutions occuring in nature are limits 
of viscous flows. It should be possible then to determine some relation- 
shio between weak solutions and solutions to the viscous flow Droblem. 
THEOREM: 

Consider the nonlinear parabolic system 

U t+ F x + B ■ AU XX (4.6) 

with initial vector U(x,0) = $(x). Here AUxx corresponds to a viscosity 
factor. Given that the strong limit (as A-+0) of the net of solution 
functions U (x,t) exists and is equal to U(x,t), then U(x,t) is a weak 

A 

solution of eqn. (4.1 ) . 

PROOF: 

Consider an arbitrary twice differentiable test vector W. Multiplying 
eqn. (4.6) by W and integrating by parts yields 



3 / 



oo 



oo 



oo 



/ / {w U A + W F(U X ) - WB) dxdt + / W(x,0)$(x)dx 

Q -oo " -oo 

t - (4.7) 

= - A / /HU, dxdt 

0 -oo x A 

Now keepina $ and W fixed, and letting A+0, the left side of eqn. (4.7) 
aoproaches the left side of eqn. (4.2), and the riaht side of eqn. (4.7) 
approaches zero. Hence U(x,t) satisfies eqn. (4.2) for all twice 
differentiable test functions. Note that U had to be a strong limit 
of (i.e. // |U^ - U| -* 0 over any bounded region in the x,t-plane) 

in order that F(U x )-*-F(U). If U-^-HJ only in the weak sense (i.e. 

U x (x,t) U(x,t) v (x,t))there is no quarantee that F(U^) would converge 
to F(U). Now it must be shown that eqn. (4.2) is satisfied for any 
arbitrary test vector W which is continuously once differentiable. Lax 
states that since U satisfies eqn. (4.2) for all twice differentiable 
test vectors W, a fortiore, U satisfies eqn. (4.2) for all once different- 
iable test vectors.^ This author disagrees with this argument since the 
class of once differentiable functions is laraer than the class of twice 
differentiable functions. Though possibly true, the statement requires 
proof. Several approaches were attempted and proved to be unsuccessful. 
One apparent way to eliminate the difficulty is to restrict our attention 
to just twice differentiable test vectors in our original definition of a 
weak solution. 

Having shown for a particular case, that the limit of viscous flow 
(as A-*0) is a weak solution. Lax conjectures that all viscosity methods 



Lax, P. D., "Weak Solutions of Nonlinear Hyperbolic Equations and 
Their Numerical Computation," Communications on Pure and Applied 
Mathematics , v. 7, n. 163, 1954. 
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should converge to the same weak solution. However he proposes a 
different limiting process, namely a special finite difference scheme, 
which is independent of the viscosity of the fluid. This approach has 
a significant advantage over the viscosity approach in that it accurately 
represents the actual behavior of the fluid in every region of the flow. 
Lax's method replaces the differentiations by finite-difference operations 
accordina to the scheme 



f x (x,th [ f ( x+Ax » t ) " f(x-Ax,t)] 
f t (x,t)- 2Zt [f( x >t+At) - f(X+AX>t ) \ 



(4.8) 



Finally, substantial numerical evidence supports Lax's conjecture 
that when the above finite difference scheme is applied to any single 
homogeneous first order conservation law 

u t + [f(u)] x = 0 f" < 0 (4.9) 

with U(x,0) = 4 >(x) , the solution U(x,t,Ax) approaches the same weak 
solution generated by the viscosity methods. The solution is given by 

u(x,t) = g (-p 5 -) (4.10) 

where y = y (x,t) maximizes 
J o o 

/ <f,(s)ds + tG(^) (4.11) 

0 r 

and 

f'Cg(s)] = s ; G'(s) = g (4.12) 

In most physical situations U(x,t) is a piecewise differentiable function. 
Then it is an easy matter to show U(x,t) is a weak solution to ean. (4.9) 
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since verification of eqn. (4.2) may be avoided. All that must be shown 
is that U(x,t) satisfies eqn. (4.9), wherever it has well-defined first 
derivatives. Since the first derivatives of U(x,t) are undefined at most 
on a set of measure zero, and U(x,t) satisfies eqn. (4.9) everywhere 
else, U(x,t) must by necessity satisfy the intearal relation (4.2) 
everywhere. Hence (4.10) is a weak solution of eqn. (4.9). The inteqral 
eqn. (4.11) is necessary to uniquely define the weak solution. 
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V. CONCLUSIONS 



The original intention of this paper was to analyze difficulties in 
a numerical system called "PUFF" [Ref. 2], "PUFF" is an attemot to 
numerically compute the reaction of a multi-layered medium to violent 
shocks. The Puff code employs Von Neumann and Richtmyer's artificial 
viscosity term. It has been shown that there were very large discrep- 
ancies between the "PUFF" solution and the classical solution inside 
the shock regions [Ref. 1]. The reason for these errors should now be 
apparent. Von Neumann and Richtmyer's artificial viscosity term had 
to comply with certain requirements. However these requirements were 
to a large degree independent of actual physical considerations and, 
hence, were not sufficiently restrictive. Consequently, in shock 
regions where viscosity is significant, the artiftctal term may not 
accurately describe the actual fluid flow. This difficulty becomes 
critical in a multi-layered medium since the number of shock waves is 
increased when the original shock reflects from several boundaries. 
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FIGURE 1 [Ref. 7] 
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